function plotrocks(filenumber)
S=readdata(0,2);
refsurface=readdata(0,25); refsurface=refsurface(:,4);
surfacehist=zeros(S,length(filenumber)); timehist=zeros(1,length(filenumber));
timeindex=1;
for count=filenumber
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
%    mx=readdata(count,10);
%    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
%    mtype=readdata(count,14);
%    mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
%     P=readdata(count,23);
Nsurfaces=readdata(count,24);
surfaces=readdata(count,25);

count
x=zeros(1,S);
for s=1:S-1 x(s+1)=x(s)+dx(s); end

sealevel=0.5*(surfaces(1,Nsurfaces)+surfaces(S,Nsurfaces));
%surfacehist(:,timeindex)=surfaces(:,4)-40e3
%surfacehist(:,timeindex)=surfaces(:,4)-refsurface;
tempsurface=surfaces(:,4)-refsurface; surfacehist(:,timeindex)=tempsurface-mean(tempsurface);  %correction for long term drift
timehist(timeindex)=time/(60*60*24*365*1e6);
timeindex=timeindex+1;
end
save surface_time timehist surfacehist x;
[TIMEHIST X]=meshgrid(timehist,x);
figure(1); clf; surf(X,TIMEHIST,-surfacehist); xlabel('Distance/m'); ylabel('Time/Ma'); set(gca,'zdir','rev');
figure(2); clf; pcolor(X,TIMEHIST,-surfacehist); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar;
end